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In this work we first study the quantum diffusion in a volume of a crystalline solid at high 
interstitial concentrations when the effects of the short-range interactions between the diffusing 
particles are to be the factors. Within the scope of the small-polaron formalism the transition rate 
depending on both the temperature and the interstitial concentration has been calculated. Then, 
using the obtained result, we consider the model of surface diffusion that reproduces qualitatively 
the diffusion behaviour of the hydrogen isotopes on W(l 10) surface, observed in the lowtemperature 
experiments [S.C. Wang and R. Gomer, J. Chem. Phys. 83, 4193 (1985)]. The coverage dependence 
of the diffusion coefficient is determined by Monte Carlo simulation. This model allows to suppose 
that substantially the different diffusion of hydrogen (tritium) and deuterium can be the outcome 
of their different quantum statistic and the direct interactions between the adparticles. 

PACS numbers:66.30.Jt, 68.35. Fx 



I. INTRODUCTION 

Quantum diffusion both in volume of solids and on solid surface has been attracting considerable attention for 
the past several decades mainly due to extraordinary high mobility of interstitial hydrogen and its isotopes at low 
temperatures. Theoretical consideration of this problem rests on low concentration of interstitial particles. This 
allows to reduce duffusion processes to tunneling dynamics of syngle particle, moving in an effective potential and 
interacting with the solid environment, in particular, with phonon sybsystem of the perfect lattice Q-|^,|^|2^]. fn 
different approaches interection with environment is assumed to be either linear or bilinear with respect to phonon 
variables and considered within the scope of perturbation theory. As a result of tunneling a particle transfers from 
one interstitial site to an adjacent that. Near T = corresponding motion is due to band- like propogation (coherent 
diffusion). At increased temperatures localization of a particle in a site occurs by interaction with phonons and the 
hopping regime is realized (incoherent diffusion). In a metal at very low temperatures the tunneling dynamics of 
particles proves to be strongly depending on nonadiabatic interection with the conduction electrons [p3 E6| . 

To calculate the diffusion coefficient different approahes have been suggested to a certain extent using the above 
physical pictures. The review of some theoretical methods and subsequent generalizations are given in the papers 



At the finite impurity concentrations the interaction between diffusing particles can distort the low concentration 
picture. So, Kagan and Maksimov have shown that when particle concentration exceeds some low critical value 
elastic interaction between particles at T — results in suppression of the band motion and in complete particles 
localization. 

At high interstitial concentrations the diffusion process is to be determined to a great extent by short-range interac- 
tions between diffusing particles || . An interesting example in this sense is the unusial diffusion of hydrogen isotopes 
on the (lfO) surface of tungsten at low temperatures revealed by DiFoggio and Gomer |^9| and Wang and Gomer 
p0| . They have found that dependence of the diffusion coefficient D on coverage 6 in all cases displays oscillator-like 
behaviour. However while D for hydrogen (and tritium) decreases at — > I , the diffusion coefficient for deuterium 
increases. This phenomena has not been explained up to now. Muttalib and Sethna |^] and Gomer ]3l) have supposed 
.£h ! that it may be a consequence of the difference in quantum statistics of 1 H and 2 H. From our point of view a double 
effect takes place here: the existence of interaction between diffusing adatoms and difference in statistics. 

In this work we first study quantum diffusion in a valume of a crystalline solids at high concentrations of interstitials, 
when the hopping regime is realized. Because of the large distortion of lattice the corresponding temperature region 
can extend to lowerest temperatures. We confine ourselves to a temperature region in which electronic degrees of 
freedom can be eleminated from consideration within of adiabatic approximation. Allowing for short-range interaction 
between interstitials, we find an expression for the transition rate depending on concentration. Then a model of surface 
diffusion , showing the same behaviour as it does for hydrogen and deuterium on VF(IIO), is considered. The Monte 
Carlo method has been utilized to simulate surface diffusion of adatoms and to calculate the coverage dependence of 
the tracer diffusion coefficient. In conclusion on the basis of the considered model we discuss the reasons for difference 
in diffusion of 1 H and 2 H atoms on VF(IIO). 
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II. TRANSITION RATE 



When considering the diffusion process in a crystal at a low interstitial concentration one proceeds from the suppo- 
sition that the jumps are sufficiently rare events, so in the intervals between the jumps a thermodynamic equilibrium 
is reached and the system vibrates around a stable crystalline configuration. The same situation occurs at high 
concentrations. In this case the mobility of particles will be determined by their jumps over the vacant interstices, 
whose concentration is small. So, between the jumps in such a crystal thermodynamic equilibrium is also reached and 
a system will vibrate around the corresponding crystalline configuration. 

Let the vibrational states of the system (host atoms plus interstitials) be described by the wavefunctions > 
depending on the configuration of the occupied interstices r — (r^r^, . . .), with {n} denoting the set of quantum 

(r) . . . (r) 

numbers. W n } take no account of the jumps and are the eigenfunctions of the some Hamiltonian Hq — T + Uq . 

For different configurations r and r , W n i and wi n >\ are localized in different regions of the configuration space and 

their overlap is small W . Considering these states as approximately ortogonal we can write down the transition rate 
as §§ 

W = < W (rr "> > r , (1) 

W (rr>) = ( 2 iT,\< O - *oK$ > f<X<? } - E$)) {n} . (2) 
{»'} 

Here and further r and r' denote the initial and final interstitial configurations during the jump of an interstitial to 
the nearest adjacent site, H = T + U is the exact Hamiltonian allowing for the jumps; < ... >{ n }i < •■■ > r are the 

thermal average over the initial states {n} and the average over all possible initial configurations r, respectivly; 

(r) 

is the energy of the state kI. 

To calculate W let us introduce the following approximations. 

(i) Taking into account that the overlap between an< ^ ^{n'} * s smau ^ we can simplify the matrix element in Eq. (^|) , 
using approximation 

< O - Ho\4J } >= A*<'> < <? } | > . (3 ) 

where A$( r ) = U— is calculated at a point of the overlap range, corresponding to some "transition" configuration 
of interstitials r* . In this configuration the jumping interstitial is situated in the middle of the way between the initial 
and the final sites during the jump. The positions of other particles for present remain not determined, ft should be 

(r) 

noted that, actually, the difference U — Uq involves only degrees of freedom related to the jumping interstitial and 
its nearest neighbours and is notably distinguished from zero in some finite region of the crystal. 

(ii) By occupying some interstice an interstitial interacts with the neighbouring host atoms. The additional force 
repulsion from the nearest interstitials must result in the sizable localization in sites even of such light particles 

(r) . . (r) 

as hydrogen. In this case one can approximate the states wi^} the eigenfunctions Hq , taken in the harmonic 
approximation 

H { o ] - < r) + \ E m k x\ + \ £ *i\ Xk - x$){xi - *£>) • (4) 

k kl 

(r) 

Here x±, Xi, ■ ■ ■ x^n are the coordinates of all N particles of the system, x k0 are the equilibrium positions of particles 

(r) (r) 

for the configuration r; <5jy are the force constants and <1>q is the potential energy of the system at the equilibrium 
position; = m is the host atom mass, if Xk relates to the lattice coordinates and = M is the interstitial mass, 
if Xk relates to the interstitial degrees of freedom. 

It is known (33) that when the volume of the crystal V — > oo, any additive quantity V - 1 A^ r \ due to the spatial 
homogeneity of the system, has the determinate limit Aq not depending on r, i.e. 

lim V~ l A^ = lim < V^A^ > r = A . (5) 

V — >OG V — >OG 

As $g r ' ) has the property of " self- averaging" (J|), taking a sufficiently large volume V, one will has to any degree of 
accuracy 
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Further we shall take $0 = 0. 

Proceeding to the normal-mode coordinates q£ , one writes down Hamiltonian (Q) as 

Here we have also allowed for the density of the vibrational states (distribution function of the frequencies of the 
normal modes) v^ r \us) has the property (Q), i.e 

lim i/W(w)= lim < u {r \uj) > r = u(ui,c), (7) 

V — >oo V — >oo 

where c is the interstitial concentration (the number of the interstitials divided by the number of the interstices). 
Thus, for a large enough crystal the frequency spectrum does not practically depend on the configuration of the 
occupied interstices r. 

In the paper |3^| it has been demonstrated that the overlap integral between the eigenf unctions of Hamiltonian (|^) , 
considered in the configuration space, one can represent as 

< ^? } >= II / d «Wn k (q k + bt' ] )Vn> k (ft), (8) 

fe J 

where <p n (qk) are the normalized harmonic oscillator eigenfunctions with the frequency Uk and the mass m and 

\b { "' ] \ = a(M/ m y\a { ^ \ ~ a(M/m)' a(M/m) I V3N . 

Here (M /m)^ a^? is the matrix connecting the displacement x s — x'fj with the normal coordinates g^; a is the jumping 
distance, a(M/m) is some dimensionless function and the interstitial is assumed to hop along the x s direction, i.e. 

I M _ (r')| _ 
FsO x s0 I — a - 

The matrix elements in (ph are proportianal to the factor J = e~ A , where A ~ (a 2 / < >o) 1 and < Ui >q is 
the mean square of the displacements of the zero vibrations interstitials. 

Allowing for approximation (||) and inserting (||) into (||) one can calculate W^- rr ' employing the small-polaron 
formalism |54|. Taking limits N — > oo , V — > oo and allowing for the fact that A<I>( r ' does not depend on either N or 
V, one obtains p5[ | 

oo 

wirrl) = \A^l e _ ST[c) r dt[eG{t . c) _ ^ 
h J 

— oo 

where 

a 2 a 2 M f , , Tiu> 



„ , . a z a z M f , . , 
dt(c) = — ^ — / aw v(lu,c) lu coth 



. . a a M f . . ?ia; 
G(i;c) = — — — / aidv(u!,c)LJ esc— — — cos Lot. 
2n J 2kbT 

Now the averaging over the different configurations r in Eq.(|]) is to be carried out in the finite region where A<£ >( - r ) 
is notably distinguished from zero. To carry out the averaging, one notes that in the first approximation A$^ r - ) is 
an anharmonic contribution in the potential energy of an interaction between the jumping particle, localized in the 

(r) 

middle of the way between the initial and the final sites during the jump, and its nearest environment. Let n k be 
the number of /c-th nearest interstitials of this particles, corresponding configuration r. Then one can write down 

A$W = AH^4 r Vf, (10) 
fe 
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where A$ is the contribution of the host atoms to A<I>( r ' and is the contribution of two-body interactions between 
the jumping particle and one of its A:-th neighbours among of the intertstitials. 

The Eq(|l(]) shows that A«3>( r ) depends only on occupation numbers rather than on placing of interstitials over 
equivalent interstices. It is easy to see that the averaging |A$( r )| 2 over the different configurations is equivalent to 

(r) 

averaging of the numbers n k over the binomial distribution with the probability p = c that the giving interstice is 

(r) 

occupied and the probability of the contrary event q = 1 — c, independly for every k. The random variables n k take 

the following values n k = 0, 1, 2 . . . n&, where nk is the maximum possible number of k-th neighbours among of the 
interstitials. 

In addition let there exist also three-body interactions between interstitials, for example, with the participation of 
the first two neighbours of the jumping particle. If is the average contribution to A$( r ') from each one of such 
three-body configurations the one can write 

A$« = A$ + J2 n k )v k 2) + -n ( {\n[ r) - l)vf 3) +..., (11) 
fe 

where dots denote the contributionss from thre-body interactions with the particuipation of the following neighbours 
and also the possible contributions from n-body interactions with n > 3. 

The Eq.(|ll|) depends on the occupation numbers only and its square has the following form 

|A$«| 2 = |A$| 2 l + ^^+^^n^^^^^^nM + J , (12) 

\ i ij ijk J 

where bi, bg, bj^, ■ ■ ■ are the known dimensionless coefficients. 

The Eq.rtlJ) can be averaged over the configurations as it has been indicated above. Taking into account that the 
average < rr >, where n are the occupation numbers of the interstices, over the binomial distribution is the polinomial 
in p = c of degree I, one obtains the transition rate W —< W^ rr ' > r 

W{T, c) = Q2n-2(c)W (T, c) , (13) 

where 

oo 

Wq(T, c) = e- s ^ f dt[e G ^ - 1] , 

h J 



Qin-i(c) = 1 + cue + a 2 c 2 + ... + a 2n -2C 2n 2 . (14) 
Index n in Eq.(H4) indicates that the interstitial's potential energy contains n-body potentials as well. 



III. MODEL OF QUANTUM DIFFUSION ON SOLID SURFACE 



The expression for transition rate obtained in the preceding section can be applyed to the study of the diffusion 
on solid surface without any restrictions. In this section we shall consider a model system , in which depending on 
parameters of the model, surface diffusion of the particles exhibits the same behaviour as the diffusion of hydrogen 
and deiterium on W(110) [^,^0|. First we calculate the frequency spectrum of the defective crystal to find W from 
Eq. (|l3"|) , employing Dean's direct numerical method |56 37 j . Then the coverage dependence of the tracer diffusion 
coefficient is determined by Monte Carlo simulations on a fixed two-dimensional lattice. 



A. Frequency spectrum 

We use a model in which the solid surface is represented by the cubic lattice of size n x I x 3, consisting of three 
layers, with regular sites occupied by substrate atoms of mass m. The adsorbed particles of mass M and with 
desired coverage 9 are localized in the interstices of the crystal lattice ranging in the upper layer. In our model the 
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particles are connected to their nearest neighbours by central and noncentral harmonic springs. This choice affords 
the simplification that the motions in each cartesian direction are independent [|36|]37[ . 

In this section we denote the displacement of the i-th particle in x-direction as Xi . The enumeration of the particles 
as well as the force constants connecting the nearest neighbours are shown in Figs. ^ and ^, where the notation is 
introduced 



[s] = ni + n 2 + ■ ■ ■ + n s -i + 2(s — l)n , 
[1] =0. 



1 < s < 1 + 1, 



The difference p s — n s — n is the number of the interstitials localized in the interstices of the s-th vertical column 
(along y-direction) . 

Let kij be the harmonic force constant connecting the nearest neghbours i and j. Then the equation of motion for 
i-th particle can be written as 



(15) 



The prime at the sum denotes the omission of the term with j = i . In addition the periodic boundary conditions 
were imposed in y-direction and rigid wall boundary conditions in the ^-direction. For simplicity it has been assumed 
that the adparticles and substrate atoms of the lower layer do not interact. Thus there exist the lateral interactions 
for adparticles only. 

_Allowing for Xi — —uj 2 Xi and introducing the notations Ui — rryJ 2 Xi, uj 2 — lo 2 /lo 2 , lo 2 — Ak/m, one has instead of 



(0) 



UJ 2 Ui 



= <\,ll, + ^2 Pij U 3 



(16) 



where 



rrii^—^ k ' 



771 kij 

Pa = — 



4k y?mim,j 



The system of the equations (16) is equalent to the eigenvalue problem 

M« u = uJ 2 u 



(17) 



for the symmetric matrix M.^ of order N x TV, depending on the configuration of the occupied interstices and having 
a block tridiagonal form 



/ Ax B 2 

a 2 b 3 



\ 



Bo A3 B4 



V 



(18) 



Bf A J 



Here B^" is the transpose of B s . The explicit forms of submatrices A s and B s are given in the Appendix. 

To calculate the spectrum of eigenvalues of the matrix M^ r ^ we use Dean's direct numerical method. Originally 
it was applied for finding the frequency spectrum of the disorded lattices with substitution impurities and has been 
discussed in detail in Refs. |56],[37|]. Here only the main poins are outlined. 

First one finds the integrated frequency spectrum N(uJ 2 ) which gives the number of the squared frequencies less 
than some uJ 2 . If ??(X) denotes the number of negative eigenvalues of a matrix X, then N(ZJ 2 ) = rj(M. < - r ' > — tU 2 I) and 
I is the unit matrix. Calculation of r\ (X) , where X is the matrix of form ( p"8| ) is based on the negative eigenvalue 
theorem by Dean and Martin p3 |. In the next step the squared frequency distribution function G(tJ 2 ) is determinated 
by 



G(td 2 



N(lu 2 +dQ 2 )- N(Q 2 ) 
du 2 N 
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After G((D 2 ) has been found the frequency spectrum v(uo) — 2luG(oj 2 )/u 2 . is calculated. The system's sizes are picked 
so that v(uj) depends only on a coverage 9 rather than the configuration of the occupied interstices. 

To calculate the frequency spectrum we consider the lattice containing 30 x 30 x 3 regular sites. The adparticles are 
distributed at random on interstices with the given probability 9. At the picked lattice's size the frequency spectrum 
practically does not depend on particular adparticle configuration. The following ratios between the force constants 
have been taken: ki/k = kijk = 1, k^jk = 8, k^/k = 9. The mass ratio for a substrate atom and adatom has been 
taken the same as for system 1 H/W, 2 H/W. The calculations have been carried out at different coverages from the 
region 0.4 < 9 < 0.9 and showed the insignificant variation of the distribution function G(uJ 2 ) at the increase of 9 in 
both cases. The typical histograms for G(tJ 2 ) are showed in Figs. || and [|. 



B. Monte Carlo simulation of surface diffusion 



For simulation of the diffusion process on the surface we consider lattice gas of the N adsorbed particles occupying 
at the regular sites of L x L square lattice. In our model Hamiltonian of the system has been taken in the form 
allowing for two and three-body interactions between adatoms 

H = Ji ^ c t c j +^2 ^ (kCj + J t ^2 CiCjC k . (19) 

(NN) (NNN) (t) 

Here J\ > and Ji < are the nearest (NN) and next - nearest neighbour (NNN) interaction energies, J t > is 
three - body interaction energy (t) . Corresponding adatom configurations are shown in Fig. |5|. The positive constant 
J n > corresponds to repulsion. Cj is the occupation number of site i and takes the value 1 and for full and empty 
site respectively. Double occupancy of sites is excluded. 

Within the scope of the random walk theory the tracer diffusion coefficient for the twodimensional system is defined 

as 

1 N 

D* = lim — V < Ar 2 (t) > , (20) 
t^oo 4Nt ^ 

i=l 

where t is the elapsed time and < Ar 2 (t) > is the mean-square displacement of the i-th particle. 

To calculate D* Monte Carlo simulation of surface diffusion has been carried out in the canonical ensemble applying 
the Metropolis importance sampling [pq| . The actual procedure was the following. An initial adatom configuration 
corresponding to the desired adatom coverage 9 was given by the random pick of 9L 2 sites. Then the temperature T 
was established and employing the standard procedure a large number of Monte Carlo steps (MCS) were performed 
to reach the thermodynamic equilibrium. The equilibrium was assumed to be reached when the total energy of 
the system started to fluctuate about its average value. Typically in a lowtemperature region about 10 3 MCS's are 
required to establish equilibrium in the lattice containing 50 x 50 sites at 0.4 < 9 < 0.9 for J\ = 2 kcal/mol, J2 = — 1.5 
kcal/mol, Jt — 1 kcal/mol . 

After thermodynamic equilibrium had been established a diffusion run was started. Here we follow the algoritm 
employed in Refs. |59],[lO]]. The elementary step of surface diffusion is assumed to be a jump of an adsorbed particle 
from the occupied initial site to an empty nearest neighbour one. First of all, the initial site of lattice is picked at 
random, if filled, an adjacent final site is randomly choosen. If this site is empty the jump can occur with some 
probability Pj, otherwise no jump occurs and a new initial site is picked. One MCS corresponds to L 2 the random 
picks of the initial site for the jump. 

The time scale is defined by the span At in which an adatom is allowed to attempt a jump to an adjacent site. If 
the mean adatom's lifetime in a site is r = 1/W , where W is given by Eq.@, the probability of a jump during At 
one can write as 



At 

T 



In order to reduce computational time we have taken At = 1/W max with W max = maxVF(T, 9) in the region of values 
T and 9 we are interested in. Thus the jump probability is given by 

Pj = ^LR. (21) 
W 

* v max 

The choice of (Ell) means that when W = W max the jump to the nearest vacant site occurs with the unity probability. 



G 



To calculate W = Q2u-2Wq we have used the frequency spectrum obtained in the preceding section and have taken 
n = 3 that corresponds to the existence of pairwise and triplet interactions of adatoms. To reproduce the behaviour 
of hydrogen isotopes in the diffusion experiments we consider two model expressions for Q4 

Q ( P{6) = 1 -0.6 9 + 6.1 9 2 - 11.8 9 3 + 5.7 9 4 (22) 

that corresponds to 1 H and 

Q { 4 2) (9) = 1 - 9.9 9 + 62.2 9 2 - 122.1 9 3 + 73.6 9 4 (23) 

for 2 H. Their behaviour is shown on Fig. |^. The shape of the curves are simular but while the minimum of is 
located in the range 6 < 1 the minimum of is displaced into the range 9 > 1 . The coefficients of the polynomials 
are chosen in such way that the functions Q\ should have the given values at the points of extremum. The remaining 
indeterminacy in the coefficiets practically does not affect the behaviour of D* . 

To calculate the diffusion coefficient D* diffusion runs of 1000 MCS's for 40 different initial adatom configurations 
were performed. The results of the calculation of coverage dependence of D* at 40K and a = 1 (see Eq.(^)) are 
presented on Fig. fj]. Figs. || and § show the temperature dependence of D* for the different coverages. 



IV. DISCUSSION 



In this work we have considered quantum diffusion in a crystalline solid at high interstitial concentrations. Our 
main result is the equation ( |l3| ) for the transition rate taking into account both the temperature and the concentration 
dependences. Futher we have studied the model system imitating quantum diffusion of adatoms on the solid surface. 
The diffusion behaviour found in this model can be used for qualitative explanation of the experimental situation 
revealed in diffusion of hydrogen isotopes on W(110) surface J29 , ||(| , 

It is known that the value measured in an experiment is the chemical diffusion coefficient D, given by Kubo-Green 
formula. For qualitative discussion we can make use the approximate expression ]3l[ | 

D <N> D* 

<(5N) 2 > ' 

where < (SN) 2 > / < N > is the normalized mean-square fluctuations of the adatom's number. In the region of 
low temperatures and high coverages the thermodynamic factor < N > / < (5N) 2 > is insignificantly altered, so the 
behaviour of D is mainly determined by the tracer diffusion coefficient D* whose coverage dependence is shown on 
Fig. uj Qualitatively the curves presented on this figure show the same behaviour as the experimental curves for 1 H 
and from Ref. @. 

The coverage dependence of D* is mainly determined by the polynomial Qa{9). The choice of Q4 = corresponds 
to 1 H diffusion and Q4 = Q4 to 2 H diffusion. As the coefficients of Q4 =< |A$' r - ) | 2 > r depend on "transition" 
configuration r* , the difference between Q4 = Q4 and Q4 = indicates that these configurations for 1 H and 2 H 
are to be different. Indeed, the configuration r* for the present is not exactly defined. The additional condition can 
consist in giving of the distance R between the neighbouring adatoms involved in the diffusion jump. It is natural 
to take R equal to the average distance between the adatoms in the states discribing particle's transition to the 
neighbouring site. The corresponding wavefunctions are to allow for a symmetry in the permutations of particles. 
Owing to the difference in the statistics the close approaches are more likely for bosons (deiterium) than for fermions 
(hydrogen and tritium). So during the jump the average distance between -"-H atoms turns out more than for 2 H. As 
a result the the configurations r* shall be different in these two cases. 

Further, the coefficients of Q4 are determined by the potentials of interactions between the diffusing particles. 
When the distance between particles is small the interaction's potentials can take very large value. Then it follows 
from the above that the coefficients of for deiterium are to be more (in magnitude) than the coefficients of 

for hydrogen. As a result polynomial is rapidly altered and starts to rise at 9 — > 1 as it is shown on FigJ^. The 
diffusion coefficient for deiterium demonstrates the same behaviour. 

Fynally we can imagine diffusion motion of 1 H and 2 H as follows. Approaching closely 2 H atoms suffer the pro- 
nounced repulsions that results in the acceleration of the diffusion process. The decisive contribution to the observed 
diffusion behaviour at the high coverages comes from three (or more)-body interactions, whose potential , is probably, 
rapidly vanishing with increasing distances between particles. So mainly 2 H atoms acquire the pronounced energies 
in the three-body interactions. It is the existence of three-body interactions that dominates the rise of D for 2 H at 
9 — > 1. At the diminished coverages the diffusion motion is mainly dominated by the two-body interaction between 
particles and results in the same behaviour of D for 1 H and 2 H as it has been revealed in the experiment. 
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APPENDIX: 



The submatrices A s and B s from Eq.dlS) have the block form 



A., 



B s 



a; 


c s 









c 





C T 


G 2 








c 


B 


:) 









(Al) 



(A2) 



The form of the submatrices from (Al) and (|A2|) strongly depends on the configuration of the occupied interstices, 
so here we shall show their structure only. 

A s is a square matrix of five diagonal banded form with two additional elements in the upper right and lower left 
corners 



/ "[s] + l P[s]+2 7[s]+3 



a: = 



p \ 



v P 



l[s]+n s P[s\+n a a [s]+n s / 



where f3 = ~kx/4k Pi = 7» = A-2,i- 

C s is a matrix of order n s x n. Its nonzero elements are on p s diagonals 



/ r (s 



(s) 



\ 



c (s) 



V o 



c {s) 



where c\f = P [s]+h[s]+ris+r 

The matrix Gi (i = 1, 2) is of tridiagonal form, of order nx n with additional elements in the upper right and lower 
left corners. 



G, 



P on p 



\p 



p \ 



■ p 

P cn J 



If the index s does not relate to the side i.e. (s ^1,0 , then 

fix = 7^:(fci + k 2 + k), 



1 A)2 

= 2l (fcl + Y + fc) ' 



On sides oti are to be altered insignificantly. 
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The matrix C = — (fc2/4fc)I is of order n X n . The matrix B^, is of order n s x n s +i. Its nonzero elements are on 
p s + p s+1 diagonals. 



( & 



/? (s) 



8 {s) 



\ 



/3 (s) 



where = 

Finally B = (1/4)1 is the matrix of order 



n x n. 
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FIG. 1. Enumeration of the particles on the crystalline surface. The harmonic force constants connecting the nearest neigh- 
bours are showed by arrows. 

FIG. 2. Enumeration of the particles in the layers perpendicular to the surface. &2 is the force constant connecting the 
substrate atoms of upper and lover layers. 



FIG. 3. Squared frequency spectrum for 30 x 30 x 3 lattice at 8 — 0.7, corresponding to 1 H/W system. 



FIG. 4. Squared frequency spectrum for 30 x 30 x 3 lattice at 8 — 0.7, corresponding to 2 H/W system. 



FIG. 5. Schematic representation of the interaction constants. 



FIG. 6. Coverage dependences of Q\ and Q\ ' . 



FIG. 7. Coverage dependence of the tracer diffusion coefficient D* at 40K for Ji = 2 kcal/mol, J2 = —1, 5 kcal/mol, J t = 1 
kcal/mol; □ corresponds to 1 H/W system, O to 2 H/W system. 

FIG. 8. The temperature dependence of the tracer diffusion coefficient D* , corresponding to 1 H/W system, at the different 
coverages. 

FIG. 9. The temperature dependence of the tracer diffusion coefficient D* , corresponding to 2 H/W system, at the different 
coverages. 
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